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Abstract. After briefly describing the anelastic approximation and Glatz- 
maier's code, we present results from our preliminary studies of core con- 
vection during the hydrogen burning phase of a 15M star, as well as our 
most recent results concerning convection in the oxygen shell of a 25M 
star. 

1. Introduction 

Two-dimensional studies of turbulent convection in the interiors of massive stars 
have been quite common (Deupree 2000; Asida & Arnett 2000; Bazan &: Arnett 
1998), but the problem is inherently three-dimensional and many aspects of 
the convective flow will not be reproduced correctly in two-dimensional models 
(Canuto 2000). Only recently has it become possible to begin tackling the 
problem of simulating convection in all its three-dimensional g(l)ory, and several 
major projects are under way to do just this (Dearborn et al. 2001; Brun 2002; 
Porter & Woodward 2000). The price of this additional dimension is additional 
computational complexity, the need for faster computers and more data storage 
capabilities, not to mention difficulties in visualizing and analyzing the resulting 
data. In this paper we describe three-dimensional simulations based on anelastic 
hydrodynamics, which we hope capture the essential phenomena of the problems 
at hand, while allowing more rapid progress at a lower computational cost than 
comparable studies using fully compressible codes. 

This paper is divided into three parts. After very briefly describing the 
anelastic approximation and Glatzmaier's implementation of it in Section 2, we 
present our hydrogen core convection model and results of this simulation in 
Section 3, and the oxygen shell convection model and results in Section 4. 

In addition to the static plots shown here, please look at http://www.supersci.org 
for movies of these models. 



2. The Numerical Method 

The simulations presented in this paper were performed using anelastic hydro- 
dynamics. This approximation allows a background density stratification, but 
filters out short-period, acoustic oscillations, thereby allowing larger timesteps 
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than fully compressible codes. We believe that these attributes make anelastic 
hydrodynamics an excellent method for studying stable convective regions in 
massive stars. 

2.1. The Anelastic Approximation 

For completeness we have included the anelastic equations below. Any interested 
reader can find a much more detailed exposition of the anelastic approximation 
in Glatzmaier (1984) and references therein. 

Barred symbols denote reference state quantities, all others perturbations. 

V-(pv) = (1) 

<9v P - fdp\ 

p— = -V-(pv®v)-pV(- + * g )- (^J sgr + 2pvx« 

+V- (V(E-i(V-v)T)) (2) 

— ds 

pf— = V • (KpTVs) - TV ■ (psv) + p(e nuc + enuc) [+p(e„ + e v )\ (3) 



Here E is the rate of strain and I the identity tensor. A vanishing coefficient of 
bulk viscosity is assumed for the viscous force calculation, and viscous heating 
is neglected. 

Of particular interest is the mass conservation equation in which the Eule- 
rian time derivative of the density perturbation drops out. This in effect filters 
out sound waves and greatly ameliorates the Courant condition, which in turn 
allows larger timesteps. Note that by design anelastic codes cannot model any 
shocks, or (super)sonic flow in general. The anelastic approximation attempts 
to solve for an average solution of the flow, ignoring the propagation of any small 
scale pressure fluctuations. 



2.2. Glatzmaier's Implementation 

The numerical implementation of these equations has been performed by Glatz- 
maier, and a detailed description of the code can be found in Glatzmaier (1984). 
The code is fully spectral, with an expansion in spherical harmonics to cover the 
angular variations and in Chebyshev polynomials for the radial dependence. 

The models described here have no radiative energy transport, so there 
are only two ways to transport energy: convective and diffusive heat flux. The 
nature of the diffusive heat flux is explained below in Section 2.5. 



2.3. The Reference State 

One of the most important steps in producing these simulations is the generation 
of the reference state. We have incorporated a general equation of state and 
have added temperature and density dependent energy generation and neutrino 
losses. 
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Radial Structure As the reference state is one-dimensional, we take advantage 
of our ID stellar evolution code KEPLER (Weaver, Zimmerman, & Woosley 
1978) and fit its density and pressure to a polytrope (P = Kp 1+1 / n ). For 
the radial dependence of these quantities we solve the resulting Lane-Emden 
equation. 

Equation of State The background temperature T as well as the necessary ther- 



able Helmholtz code (Timmes & Swesty 2000). 

Energy Generation and Neutrino Losses It is important to capture the full 
temperature and density dependent nature of the energy generation within the 
convective region, and in the case of the oxygen shell also of the neutrino losses. 
For this reason we have added not only radially varying reference state nuclear 
energy generation and neutrino loss rates, but have also kept the perturbations 
to allow for differential heating between cold and hot convective bubbles. The 
nuclear energy generation rate in the relevant zones is fit to the conventional 
£mic oc pT n law, and the neutrino loss rate to e v oc p m T n . 

2.4. Boundary Conditions 

The top and bottom boundaries are modeled as impermeable and stress-free 
(v± = and (J^j = 0). This means that presently we are not able to study 
important phenomena such as convective overshoot. 

Generally we fix the entropy gradient at the bottom boundary by specifying 
a base luminosity, taken from KEPLER. For the boundary condition at the top 
we either force the entropy perturbation to vanish or we specify the top luminsity, 
i.e. the entropy gradient at the outer boundary. 

2.5. Turbulent Diffusion 

The diffusion of momentum and entropy by small-scale eddies is expected to 
dominate the diffusion due to radiative and molecular processes. The effects of 
all unresolved motions are parametrized as turbulent diffusion. In this model 
we have chosen constant turbulent viscous and thermal diffusivities, which we 
have tried to minimize to the extent permitted by the finite resolution. 



3. Hydrogen Core Convection 

As our initial testcase for this code, in order to get familiar with it, we chose 
to model core convection of a 15Mq star half-way through the main sequence. 
This stage of the star's life is not very violent; there are no neutrino losses, 
the energy generation doesn't have a very steep temperature dependence, there 
are no shocks, and density and pressure contrasts are expected to be small - a 
perfect case for an anelastic model. We view this model as a proof of principle, 
and as preparation for modeling more advanced stages of stellar evolution 



modynamic derivatives 
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Figure 1. 15 Mq hydrogen core convection model reference state. 
3.1. Parameters 

The dependence of the reference state radius, density, temperature, and pres- 
sure on the Lagrangian mass coordinate are shown in Figure 1. Table 1 below 
summarizes the remaining parameters of interest. -Mb ase is the mass of the cen- 



Table 1. Parameters of the 15Mq model of core convection during 
main sequence hydrogen burning. 



M base 


M top 


^base 


^base 


^top n e,nuc 


O.OliW© 


4.16M© 


1OOOL 


0.154^0 


1.18% 15.387 



tral portion of the star that is cut out, and M^ p is the total mass within the 
convective core of this 15Mq star. Even without the central region we are still 
modeling 99.8% of the convective core, both by volume and mass. 

3.2. Results 



Table 2. Numerical results of 15Mq hydrogen core convection sim- 
ulations. 

Resolution: n max = 241, / max = 63,m m ax = 31 



n 


"max 


< v > 


< ±pv' 2 > 


^Ra 




0.0^ 


1.5^ 


0.98^P 


2.0 x 10 10 


4 x 10 7 


1700 


6.3 x l(T 5 £f^ 


3.5^ 


1.0^ 


2.5 x 10^ 


2 x 10 7 


3000 



The Non-rotating Case The non-rotating model was run with a final resolution 
of 241 Chebyshev polynomials and spherical harmonics up to I = 63, m = 31. 
The model simulates about 44 days in the life of the star. An average convective 
speed of ~ 10 5 cm/s implies a convective turnover time of r^ = 7x 10 5 s ~ 8days, 
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Figure 2. 15 Mq hydrogen core convection model - non-rotating, 
a) entropy perturbation in equatorial slice, b) entropy perturbation in merid- 
ional slice, c) radial component of velocity at constant radius surface (68%, 
R = 5.6 x 10 10 cm) in equal-area projection, 
a) and b) also show the planar velocity field 




Figure 3. 15 Mq hydrogen core convection model - rotating, 
a) latitudinal component of velocity in equatorial slice with velocity field, b) 
kinetic energy in meridional slice, c) entropy perturbation at constant radius 
surface (68%, R — 5.6 x 10 10 cm) in equal-area projection. 
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and hence about 6 convective turnovers. The first row of Table 2 summarizes the 
numerical results of this simulation, and Figure 2 provides some visualization 
of the convective motions. Note that without rotation there is no preferred 
orientation of the model, which explains why figures 2a) and 2b) look so similar. 
Both slices seem to indicate a somewhat dipolar flow field, with high entropy gas 
rising on one side and falling at the other. Such dipolar flow pattern have been 
observed in fully compressible calculations as well (Jacobs, Porter, & Woodward 
1998). 

All thermodynamic contrasts ((Sp 1 / p)max) (5P'/P) maxj etc.) are of the 
order 10~ 5 , comfortably below 1%. The mass-averaged reference state sound 
speed is about 700km/s, so the flow is completely subsonic. We conclude that 
this problem is very suited for an anelastic treatment, as expected. 

The Rotating Case The rotating model was run for about 120 days in the star's 
life, with the same resolution as the non-rotating case. This amount of star 
time corresponds to about 15 convective turnovers. As is immediately obvious 
by looking at Figure 3, rotation has a significant effect on the convective flow 
pattern. The rotational axis breaks the symmetry of the problem, so the flow 
behaves differently in the equatorial plane than in any meridional plane and the 
dipolar flow pattern is no longer visible. In the meridional slice the strong, high 
kinetic energy flows are concentrated close to the axis, with convection being 
strongly suppressed in the equatorial regions further away from the axis. We 
also see evidence of Taylor columns, elongated, columnar flow structures parallel 
to the rotation axis. 

4. Oxygen Shell Convection 

Of the advanced stages of stellar evolution the oxygen shell burning phase is of 
particular importance. For one, convection during this brief time before core 
collapse might produce an asymmetry, leading to non-uniform mixing and per- 
haps even an asymmetric supernova explosion. Several two-dimensional studies 
of this problem have been conducted (Arnett 1994; Bazan & Arnett 1998; Asida 
& Arnett 2000), and several of these simulations found close to sonic flow, large 
density and pressure contrasts, and hence significant acoustic fluxes. 

As an alternative to these simulations, we present results of our anelastic 
model of a similar problem - oxygen shell convection in a 25Mq star, shortly 
after core oxygen depletion. As shown below, we have found a solution within 
the anelastic approximation, with small density and pressure contrasts, and 
subsonic flows. 

4.1. Parameters 



Table 3. Parameters of the 25Mq model of convection during oxygen 
shell burning, about 4 days before core collapse. 
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Figure 4. 25 Mq oxygen shell convection model reference state. 

During the modeled phase the rate of nuclear energy generation exceeds the 
rate of neutrino losses by about a factor two. The net energy generation rate 
is positive, f^ c 4irr 2 p(r){e nuc (r) + e u (r)) = 1.5 x 10 45 erg/s. KEPLER indicates 
that this energy is going into PdV work and into change of the gravitational 
potential. In the models presented here we have chosen to allow this excess 
energy to escape off the top, as we are currently unable to model expansion. 

4.2. Results 

Table 4. Numerical results of 25Mq oxygen shell convection simu- 
lations. 

Resolution: n max = 145, l max = 63, m max = 31 







< v > 


< ±pv 2 > 


^Ra 


^Rc 




180^ 


49^£ 


1.7x10^ 


5 x 10 7 
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0.08^ 


170^ 


40^P 


1.2 x 10"^ 
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Table 5. Thermodynamic contrasts for the 25Mq oxygen shell con- 
vection simulations. 





(8jf\ fM\ 
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\ 1 / max 
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6 X 10" 4 
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3 x 10" a 3 x 10~ a 


1 x 10" a 


1 x 10~ 3 



The Non-rotating Case The non-rotating model was run for an equivalent star 
time of 6500 seconds. Given an average speed of about 50 km/s we determine the 
convective turnover time to be about 70 seconds, which implies that we followed 
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Figure 5. 25 Mq oxygen shell convection model - non-rotating, 
a) entropy perturbation in equatorial slice with velocity field, b) kinetic energy 
in meridional slice, c) radial component of velocity at constant radius surface 
(84%, R — 7.4 x 10 8 cm) in equal-area projection. 



■ ■ 




Figure 6. 25 Mq oxygen shell convection model - rotating, 
a) entropy perturbation in equatorial slice with velocity field, b) radial com- 
ponent of velocity in meridional slice, c) kinetic energy at constant radius 
surface (68%, R — 6.8 x 10 8 cm) in equal-area projection. 
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Figure 7. 25 Mq oxygen shell convection model. Left: Comparison 
of convective velocities with sound speed. Right: Differential rotation 
in equatorial plane of the rotating model. 



convection for about 90 turnovers. Figure 5 shows that the convective speed 
remains around two orders of magnitude below the sound speed throughout the 
entire modeled region. All thermodynamic contrasts are below the 1% level, 
as shown in Table 5. It is important to note that nothing in the model forces 
the resulting convective velocities to be subsonic and nothing prevents the ther- 
modynamic perturbations from rising to a level comparable to the background. 
The fact that this solution remained subsonic and with small thermodynamic 
contrasts implies to us that the anelastic approximation is a valid one for the 
problem of the convective oxygen shell. 

The entropy perturbations in the equatorial plane (see Figure 6) indicate a 
flow pattern dominated by four large cells. High entropy fluid is rising towards 
90° and 270°, and low entropy fluid is falling along the 0° and 180° directions. 
We have not yet determined how much of this simple structure is dependent on 
Rayleigh and Reynolds numbers. It is possible that higher resolution simula- 
tions, with higher Rayleigh and Reynolds numbers, will destroy this pattern. 

The Rotating Case The rotating case was run for 3300 seconds, which corre- 
sponds to about 50 convective turnovers. Again the convective speed is very 
subsonic and the thermodynamic contrasts are still below 1%, albeit somewhat 
higher than in the non-rotating case. 

Although the reference state is rotating with a fixed angular speed (Q = 0.08 
rad/s) the flow is free to arrange itself into any form of differential rotation. The 
result can be seen in Figure 5. Our rotating model shows a lot of differential 
rotation, with the outer regions of the equatorial plane rotating nearly twice 
as fast as the inner ones. Just like in the hydrogen core convection model we 
see Taylor-column-like structures parallel to the rotation axis in the meridional 
planes. 

5. Conclusions 

Anelastic hydrodynamics appears to be a useful tool for studying convection in 
the interiors of massive stars in three dimensions. As long as the problem at 
hand satisfies the anelastic conditions of subsonic flow speeds and small ther- 
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modynamic perturbations, a lot of progress can be made with an anelastic code 
in a shorter amount of time and with less computational resources than with a 
fully compressible code. 

In this paper we have presented results of our efforts to model convection 
in the core of 15Mq main sequence star as well as in the oxygen burning shell of 
a 25Mq star at a much more advanced stage of stellar evolution. The resulting 
convective speeds are in agreement with results from KEPLER using mixing 
length theory. The resulting three dimensional flow pattern are not always 
uniform, and can show strong asymmetries, such as the dipolar flow in the 
hydrogen core convection model. We have found a significant difference in the 
flow pattern of our rotating models compared to the non-rotating ones. The 
presence of rotation breaks the spherical symmetry and allows more vigorous 
convection along the axis of rotation. The rotating oxygen shell model showed 
a significant degree of differential rotation. 

This work was supported by the Scientific Discovery though Advanced Com- 
puting (SciDAC) Program of the DOE (DE-FC02-01ER41176). 
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